Elastic Theory of Defects in Toroidal Crystals 
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We report a comprehensive analysis of the ground state properties of axisymmetric toroidal crys- 
tals based on the elastic theory of defects on curved substrates. The ground state is analyzed as 
a function of the aspect ratio of the torus, which provides a non-local measure of the underlying 
Gaussian curvature, and the ratio of the defect core-energy to the Young modulus. Several struc- 
tural features are discussed, including a spectacular example of curvature-driven amorphization in 
the limit of the aspect ratio approaching one. The outcome of the elastic theory is then compared 
with the results of a numerical study of a system of point-like particles constrained on the surface 
of a torus and interacting via a short range potential. 



I. INTRODUCTION 

The study of the spontaneous organization of micro- 
scopic objects such as colloids, amphiphiles or protein 
clusters into mesoscale structures (mesoatoms) is an ac- 
tive research area that offers challenges on many fronts^'^. 
One would like to create a rich warehouse of raw materi- 
als from which to engineer mesomolecules or bulk materi- 
als with novel mechanical, optical or electronic behavior. 

A fertile source of potential mesoatoms is provided 
by the self-assembly of micron or smaller scale colloidal 
particles on curved interfaces. The interplay between 
spatial curvature and condensed matter order in these 
systems has proven to be very rich. Spatial curvature 
leads to novel defect arrays in the ground state that have 
widespread implications for the fundamental and applied 
physics of curved phases of matter with crystalline, hex- 
atic or nematic order. 

In a recent paper Kim et al reported the formation 
of toroidal micelles from the self-assembly of dumbbell- 
shaped amphiphilic molecules^. Molecular dumbbells 
dissolved in a selective solvent self-assemble in an aggre- 
gate structure due to their amphiphilic character. This 
process has been observed to yield coexisting spherical 
and open-ended cylindrical micelles. These structures 
change slowly over the course of a week to toroidal mi- 
celles which thus appear more stable. Toroidal geome- 
tries also occur in microbiology in the viral capsid of the 
coronavirus torovirus'^. The torovirus is an RNA viral 
package of maximal diameter between 120 and 140 nm 
and is surrounded, as other coronaviridae, by a double 
wreath/ring of cladding proteins. 

Carbon nanotori form another fascinating and tech- 
nologically promising class of toroidal crystals^ with re- 
markable magnetic and electronic properties. The inter- 
play between the ballistic motion of the tt electrons and 
the geometry of the embedding torus leads to a rich va- 
riety of quantum mechanical properties including Pauli 
paramagnetism^ and Aharonov-Bohm oscillations in the 
magnetization^. Ring closure of carbon nanotubes by 
chemical methods^ suggest that nanotubes may be more 
flexible than at first thought and provides another tech- 
nique of constructing carbon tori. 

A unified theoretical framework to describe the struc- 



ture of toroidal crystals is provided by the elastic theory 
of defects in a curved background^" "'^^. This formalism 
has the advantage of far fewer degrees of freedom than a 
direct treatment of the microscopic interactions and al- 
lows one to explore the origin of the emergent symmetry 
observed in toroidal crystals as the result of the inter- 
play between defects and geometry. The latter is one 
of the fundamental hallmarks of two-dimensional non- 
Euclidean crystals and leads to universal features ob- 
served in systems as different as viral capsids and carbon 
macromolecules . 

In this paper we provide a detailed analysis of the 
structural properties of toroidal crystals. We show that 
the ground state has excess 5— fold disclination defects 
on the exterior of the torus and 7— fold defects on the 
interior. The precise number of excess disclinations, as 
well as their arrangement, is primarily controlled by the 
aspect ratio of the torus. Since defective regions are 
physically distinguished locations they are natural places 
for biological activity and chemical functionalization. A 
thorough understanding of the surface structure of crys- 
talline assemblages could represent a significant step to- 
wards a first-principles design of entire libraries of nano 
and mesoscale components with precisely determined va- 
lence. Such mesoatoms could serve in turn as the building 
blocks for mesomolecules or bulk materials via sponta- 
neous self-assembly or controlled fabrication. 

The paper is organized as follows. In Sec. II we review 
the elastic theory of defects in two-dimensional curved ge- 
ometries and derive the ground state energy of a toroidal 
crystal. In Sec. Ill we summarize some fundamentals 
of the geometry of triangulated tori and discuss how the 
intrinsically discrete problem of crystallography can be 
reconciled with the result of the continuous elastic theory 
presented in the previous section. Sec. IV is devoted to 
the analysis of the crystalline structure arising from the 
solution of the elastic problem. In Sec. V we discuss the 
results of numerical minimization of the potential energy 
of a system of classical particles interacting via a short- 
range potential on the surface of a torus in the light of the 
results of Sec. IV. Finally, in Sec. VI we specialize our 
analysis to the case of a "fat torus" of aspect ratio one 
and we show how the curvature singularity at the center 
of the torus is responsible for a remarkable curvature- 




nations can be expressed as 



FIG. 1: (Color online) The standard parametrization of a 
circular torus of radii Ri and R2 . 



driven transition to a disordered state. Some of the re- 
sults presented here have been already announced in Ref. 
13. 



II. DEFECTS AND GEOMETRY 

A two-dimensional torus of revolution T^ is described 
in parametric form by: 



X = (i?i + i?2 COS Ip) COS (j) 

y = {Ri + i?2 cos tp) sin </> 
z = R2 sin ip 



(1) 



where i?i > R2 are the two radii of the torus. The metric 
tensor gij (with determinant g) and the Gaussian curva- 
ture K can be written respectively as: 
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Within the elastic theory of defects on curved surfaces 
the original interacting particle problem is mapped to a 
system of interacting disclination defects in a continuum 
elastic curved background. Disclinations are character- 
ized by their topological or disclination charge, qi, repre- 
senting the departure of a vertex from a the 6— fold coor- 
dination of a perfect triangular lattice. Thus g^ = 6 — q, 
where Ci is the coordination number of the ith vertex. 
A classic theorem of Euler requires the total disclination 
charge of any triangulation of a two-dimensional Rieman- 
nian manifold M to equal 6xm, where xm is the Euler 
characteristic of M. In the case of the torus xm = 0, 
and thus disclinations must appear in pairs of opposite 
disclination charge (i.e. 5— fold and 7— fold vertices with 
Qi — I and —1 respectively) in order to ensure disclina- 
tion charge neutrality. 

The total free energy of a toroidal crystal with N discli- 
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where Y is the 2D Young modulus and r(a;) is the so- 
lution of the following Poisson problem with periodic 
boundary conditions: 



^gT{x) = Ypix) , 
where Ag is the Laplace-Beltrami operator 
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and p{x) is the total topological charge density 
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of TV disclinations located at the sites Xk together with 
the screening contribution due to the Gaussian curvature 
K{x) of the embedding manifold. 

The first term in Eq. (4) represents the long-range elas- 
tic distortion due to defects and curvature. The second 
term in Eq. (4) is the defect core-energy representing 
the energy required to create a single disclination defect. 
This quantity is related to the short-distance cut-off of 
the elastic theory and is proportional to the square of the 
topological charge times a constant Ec^'*- Finally Fq is the 
free energy of a flat defect-free monolayer. The Gauss- 
Bonnet theorem-'^^ requires the total topological charge 
to be zero on the torus: 



N 

Y,<lk= / (fxK{x) = Q 
fe=i '' 



(8) 



The function r(a;) can be expressed in the Green form: 
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where GL{x^y), the Laplacian Green function on the 
torus, satisfies the equation: 
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(10) 



with Sg{x,y) the Dirac delta function on the torus and 
A = 47r^i?ii?2 is the surface area. It is easy to prove 
that the solution of the traditional Green-Laplace equa- 
tion doesn't exist on closed manifolds like the torus, from 
which the extra term A~^ appearing in Eq. (10). A"^ 
is the eigenfunction of the Laplacian associated with the 
null eigenvalue (zero mode) and is necessary because a 
pure isolated source (giving rise to the delta-function on 
the right-hand side) has no place to terminate the field 
on the closed torus. 



As usual the calculation of the Green function can be 
remarkably simplified by conformally mapping the torus 
to a domain of the Euclidean plane via a suitable system 
of isothermal coordinates. Intuitively the torus is confor- 
mally equivalent to a rectangular domain described by a 
system of Cartesian coordinates. To make this explicit, 
one can equate the metric of the torus in the coordinates 
('0, (j}) to a conformally Euclidean metric in the coordi- 
nates (C,??): 

ds'^ == Rldij? + {Ri + R2 costl)fd(fp- = w {df + drj^) , 

where w is a positive conformal factor. Taking rj = (f> and 
w — (i?i -|- i?2 cosip)"^, the coordinate ^ is determined by 
the differential equation: 



dip r + cos ip ' 



(11) 



where r = R1/R2, the aspect ratio of the torus, may 
be taken greater or equal to one without loss of general- 
ity. Choosing the plus sign and integrating both sides of 
Eq. (11) wc find: 
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Taking tp E [— 7r,7r], the integral (12) yields: 
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where 
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In the transformed coordinate system (^, rf) the modified 
Green-Laplace equation reads: 



AGLix,y) =6{x,y) - — , 



(14) 



where A and S are now the Euclidean Laplacian and delta 
function. The function Gl{x, y) can be expressed in the 
form: 

Gl{x, y) = Go(a;, y)-{G^{x, ■ ))-(Go(- , y)) + (Go(- , • )) , 

where Go{x, y) is the Laplacian Green function on a pe- 
riodic rectangle and the angular brackets stand for the 
normalized integral of the function Go (a;, y) with respect 
to the dotted variable: 



{Go{x.-)) 



A 



GQ{x,y). 



(15) 



Analogously the function (Go(- , • )) is given by 
{Go{-,-))= f^-^Goix,y) 



and ensures the neutrality property: 

Jd^xGL{x,y)= Jd^yGL{x,y) 



= 0. 



(16) 



Using standard analysis the Green function GQ{x,y) 
at the points x = {^jV) ^'^d y = {CiV') oi a periodic 
rectangle of edges pi and p2 can be calculated in the 
form: 



Goix,y) = 
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where z = S, + it], z' = ^' + irj' and i?i is the Jacobian 
theta function and reflects the double periodicity of the 
torus'^". A pedagogical derivation of the Green function 
Go(a;,y) is reported in Appendix A. Substituting Eq. 
(17) in Eq. (9) with p2 = 27r and 
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d^P 



- COSl/i 



KTT , 



we obtain: 
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where rs(a;) represents the stress field due to the Gaus- 
sian curvature of the torus and is given by: 
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The function rd{x,Xk) is the stress field at the point x 
arising from the elastic distortion due to a defect at Xk 
and is given by 
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where Li2 is the usual Eulerian dilogarithm and 



= yr^ — 1 — ' 



(20) 
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A derivation of the functions Ts{x) and Td{x, Xk) is given 
in Appendix B. Integrating the function r(a;) on the 
torus gives the elastic energy of an arbitrary collection 
of disclinations. A detailed analysis of the crystalline 
structure arising from Eqs. (4), (18), (19) and (20) is 
carried out in Sec. III. 




FIG. 2; (Color online) Construction of a defect-free triangu- 
lation of the torus. On top planar map of the triangulated 
torus corresponding to the (n, m, I) configuration (6, 3, 1). On 
the bottom (6, 3, 6) chiral torus. The edges of each one of the 
six tubular segments has been highlighted in red. 



III. GEOMETRY OF TOROIDAL DELTAHEDRA 

Before analyzing the defect distribution arising from 
the elastic energy of Eq. (4), together with Eq. (18), it 
is necessary to understand the geometry of triangulated 
tori. Reconciling the predictions of a continuum elastic 
theory with the intrinsically discrete nature of crystallog- 
raphy requires an understanding of the possible lattices 
that can be embedded on the torus and the associated 
defects. The problem of classifying the possible triangu- 
lations of the 2— torus has received considerable atten- 
tion from mathematicians, physicists and chemists over 
the past twenty years. Lavrenchenko^^ proved in 1984 
that all the triangulations of the torus can be generated 
from 21 irreducible triangulations by certain sequences of 
operations called vertex splitting*^. After the discovery 
of carbon nanotubes in 1991 and the subsequent the- 
oretical construction (later followed by the experimen- 
tal observation) of graphitic tori, many possible tessel- 
lations of the circular torus have been proposed by the 
community^^^^^. In this section we review the construc- 
tion of a defect-free triangulated torus and we show how 
the most symmetric defective triangulations can be gen- 
erally grouped into two fundamental classes correspond- 
ing to symmetry groups Dnh and Dnd respectively. 

For the sake of consistency with the existing literature 
we adopt here the language developed to describe the 
structure of carbon nanotubes. The structure of a trian- 
gulated cylinder can be specified by a pair of triangular 



lattice vectors c and t, called the chiral and translation 
vector respectively, which together define how the planar 
lattice is rolled up. In the canonical basis a^ — (1, 0) and 
a2 = (^, -^), the vector c has the form: 



c = nai + ma2 



meZ. 



(22) 



The translation vector t, on the other hand, can be ex- 
pressed as an integer multiple 



t = let 



l£ Z 



(23) 



of the shortest lattice vector e^ perpendicular to c. The 
vector Cj is readily found to be of the form: 



et 



{n + 2m)ai — (2n + m)a2 
[n + 2rn : 2n + m) 



where (a : 6) denotes the greatest common divisor of 
a and b and enforces the minimal length. The three- 
dimensional structure of the torus is obtained by con- 
necting the edge OT of the rectangle in Fig. 2 (top) to 
TTC and TJC to TTT. The edge OT is then mapped to 
the external equator of the torus while the edge OC to 
the (j) — Q meridian. The resultant toroidal lattice has 
characteristic chirality related to the initial choice of the 
vector c. In the nanotubes literature armchair referes to 
the lattice obtained by choosing n = m, zigzag to that 
obtained for m = and chiral to all other lattices. An 
example of a {n,m,l) chiral torus is shown in Fig. 2 (bot- 
tom) for the case n = 6, m = 3 and Z = 6. The chirality 
is extremely important in graphitic carbon nanotube or 
nanotori, where it determines whether the electronic be- 
havior of the system is metallic or semiconducting. 

By Euler's theorem one can prove that the number 
of triangular faces F and the number of vertices V oi a 
triangular toroidal lattice is given by: 
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Denoting An the area of the rectangle with edges c and t 
and At the area of a fundamental equilateral triangle, the 
number of vertices of a defect-free toroidal triangulation 
is then: 
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The planar construction reviewed above allows only lat- 
tices with an even number of vertices. Defect-free 
toroidal deltahedra with an odd number of vertices are 
also possible and their construction is generally achieved 
by assembling congruent octahedral building blocks. An 
example of this scheme will be briefly discussed in Sec. V 
for the case V — SI and r = 6. We refer the reader to 
Ref. 22 for an comprehensive review of the topic. 

The embedding of an equal number of pentagonal and 
heptagonal disclinations in the hexagonal network was 
first proposed by Dunlap in 1992 as a possible way to in- 
corporate positive and negative Gaussian curvature into 



the cylindrical geometry of carbon tubules^'^. According 
to the Dunlap construction the necessary curvature is in- 
corporated by the insertion of "knees" (straight cylindri- 
cal sections of the same diameter joined with a kink) in 
correspondence with each pentagon-heptagon pair arising 
from the junction of tubular segments of different chiral- 
ity (see Fig. 4). In particular, a junction between a (n, 0) 
and a (to, to,) tube can be obtained by placing a 7— fold 
disclination along the internal equator of the torus and 
a 5— fold disclination along the external equator. Since 
the radii of the two segments of a junction are different 
by construction, the values of n and to. are commonly 
chosen to minimize the ratio |c(„_o)l/|c(m,m)l — n/\/3m. 
By repeating the 5 — 7 construction periodically it is pos- 
sible to construct an infinite number of toroidal lattices 
with an even number of disclinations pairs and dihedral 
symmetry group Dnh (where 2n is the total number of 
5 — 7 pairs. Fig 5). The structure of the lattice is de- 
scribed by the alternation of two motifs with crystalline 
axes mutually rotated by 30° as a consequence of the 
connecting disclination. One of the fundamental aspects 
of Dunlap's construction is that all the disclinations are 
aligned along the two equators of the torus where the 
like-sign Gaussian curvature is maximal. As we will see 
below, this feature makes these arrangements optimal in 
releasing the elastic stress due to curvature. 

Another class of crystalline tori with dihedral antipris- 
matic symmetry Dnd was initially proposed by Itoh et 
al^^ shortly after Dunlap. Aimed at reproducing a struc- 
ture similar to the Cgo fuUerene, Itoh's original construc- 
tion implied ten disclination pairs and the point group 
D^d- In contrast to Dunlap tori, disclinations are never 
aligned along the equators in antiprismatic tori, instead 
being staggered at some angular distance 5ip from the 
equatorial plane. Hereafter we will use the symbol TAn 
to refer to toroidal deltahedra with 2n disclination pairs 
and Dnd symmetry group. 

A systematic construction of defected triangulations 
of the torus can be achieved in the context of planar 
graphs^^'^^. A topological embedding of a graph in a 




FIG. 4: (Color online) Dunlap knees obtained by joining two 
straight tubular segments with (n, 0) and (m, m) chirality. 
[Courtesy of A. A. Lucas and A. Fonseca, Facultes Universi- 
taries Notre-Dame de la Paix, Namur, Belgium]. 




FIG. 5: (Color online) Five-fold polygonal torus obtained by 
joining (5, 5) and (9, 0) tubular segments via ten pairs of 5 — 7 
rings. This structure was originally proposed by the authors 
of Ref. 27 as a possible low-strain configuration for carbon 
nanotori. [Courtesy of A. A. Lucas and A. Fonseca, Facultes 
Universitaries Notre-Dame de la Paix, Namur, Belgium]. 





FIG. 3: (Color online) Voronoi lattices of a TPn prismatic 
(top) and TAn antiprismatic (bottom) toroids with _Ri = 1 
and R2 = 0.3. 



two-dimensional manifold corresponds to a triangulation 
of the manifold if each region of the graph is bounded 
by exactly three vertices and three edges, and any two 
regions have either one common vertex or one common 
edge or no common elements of the graph. The sim- 
plest example of toroidal polyhedra with Dnd symmetry 
group, featuring only 5— fold and 7— fold vertices, can be 
constructed by repeating n times the unit cell of Fig. 7a. 
These toroidal antiprisms^^ have V = An vertices and 
can be obtained equivalently from the edge skeleton of a 
n— fold antiprism by attaching at each of the base edges 
a pentagonal pyramid and by closing the upper part of 
the polyhedron with n additional triangles. By count- 
ing the faces one finds F — 5n + 2n + n — 8n from 
which V = An. The simplest polyhedron of this family 
has V — \2 and Dy,d symmetry group (see top left of 
Fig. 6) and corresponds to the "drilled icosahedron" ob- 








FIG. 6: (Color online) First six toroidal antiprisms obtained 
by repeating the unit cell of Fig. 7. The first polyhedron on 
the left is the "drilled icosahedron" . 



tained by removing two parallel faces of an icosahedron 
and connecting the corresponding edges with the six lat- 
eral faces of an antiprism with triangular base (i.e. a 
prolate octahedron) . Starting from this family of toroidal 
antiprisms a number of associated triangulations having 
the same defect structure can be obtained by geometri- 
cal transformations such as the Goldberg inclusion^^"^^ . 
Such transformations, popularized by Caspar and Klug 
for the construction of the icosadeltahedral structure of 
spherical viruses^^, consist in partitioning each triangu- 
lar face of the original graph into smaller triangular faces 
in such a way that old vertices preserve their valence and 
new vertices have valence six. The partition is obtained 
by specifying two integer numbers {L,M) which define 
how the original vertices of each triangle are connected 
by the new edges so that the total number of vertices is 
increased by a factor T = L'^ + LM + M^. 

A general classification scheme for Dnd symmetric tori 
was provided by Berger and Avron^^ in 1995. Their 
scheme is based on the construction of unit graphs com- 
prising triangular tiles of different generations. In each 
generation, tiles are scaled in length by a factor l/\/2 
with respect to the previous generation. This rescaling 
approximates the non-uniformity of the metric of a cir- 
cular torus. 

Dunlap toroids can be obtained from unit cells such 
as those shown in Fig. 8. The geometrical properties of 
these graphs can be described in different ways. A par- 
ticularly intuitive way, in the spirit of this paper, consists 
in specifying the distances between 5— and 7— fold pairs. 
One starts by drawing the smallest convex loop passing 





(a) 



(b) 



FIG. 7: (a) Unit cell for toroidal antiprisms. 5— fold vertices 
are circled and 7— fold vertices are boxed, (b) Unit cell of 
a D„d torus in the Berger-Avron construction. The graph 
consists of four generation of tiles and the internal equator 
of the torus is mapped into the horizontal line passing to the 
mid-point between the 6th and the 7th vertex. 



through defective sites. This identifies a central polygon 
whose upper vertices {vi and V2 in Fig. 9) have degree 
five and lower vertices {v^ and vq in Fig. 9) have degree 
seven. Then calling a the distance between 5— fold ver- 
tices vi and V2, b that between 7— fold vertices v^ and 
vq and c the length of the segment V3V4 (d = a for a 
trapezoid), we can express the total number of triangles 
enclosed by the central polygon as: 

/ = 2c2 - a^ - 6^ 

Each 7— fold vertex sits at the apex of a diamond-shaped 
complex of /' = 7 triangles. Each 5— fold vertex, on 
the other hand, is at the apex of a triangular region of 
f" — {c— a+ 1)^ triangles. The graph is completed by a 
rectangle of height c — a + A and base of arbitrary length 
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FIG. 8: (Color online) Unit cells for Dunlap toroids of type 
(2, 1, 3, 1) and (3, 1, 3, 1) according to the classification scheme 
given here. Highlighted regions correspond to the central 
polygon. 





FIG. 9: (Color online) Central polygon (right) and rectan- 
gular (zig-zag) region in our construction scheme of Dunlap's 
toroids. In this example (a, b, c, d) — (2, 1, 3, 1). 



2d containing: 



/"' = 4d(c -a + A) 



triangles. The total number of vertices is 

Vg = f + f' + r + f"'l2 

= c^ -h^ + 2{c-a){c + d+l) + 8{d + 1) . (25) 

The final triangulation of the torus is obtained by re- 
peating the prismatic unit cell I times and therefore has 
V ~ IVg vertices. This scheme provides direct infor- 
mation on the arrangement of defective sites. Thus for 
instance an (a, 6,c, d) = (2,1,3,1) unit cell (see top of 
Fig. 8) has 5— fold vertices separated by two lattice spac- 
ings and 7— fold vertices by one lattice spacing. On the 
other hand the integers n and m giving the chirality of 
the two segments of the junction (n, 0)/(m, m) are given 
directly by as: 

n — c — a + A 
m — 2c — a — b . 

Thus the (2,1,3,1) cell of Fig. 8 is obtained from the 
junction between a (5, 0) and a (3, 3) tubular segment. 




(a) (b) 

FIG. 10: Unit cells for TP(2)n and TP(3)n toroids. 





FIG. 11: (Color online) (a) TP5a and (b) TP7b toroids with 
1/ = 35 and 49 obtained by repeating the unit cells of Fig. 10. 
7— fold vertices form dimers normal to the equatorial plane 
while 5— fold vertices are (a) distributed along the external 
equator or (b) form a double ring above and below the equa- 
torial plane. 



Dunlap's toroids are not the only examples of defec- 
tive triangulations of the torus with dihedral prismatic 
symmetry group Dnh. With the help of numerical sim- 
ulations (see Sec. V) we found two other classes whose 
unit cell is shown in Fig. 10. Unlike Dunlap's toroids, 
the 7— fold vertices in these prismatic triangulations are 
not aligned along the internal equator of the torus, but 
rather grouped in dimers normal to the equatorial plane. 
5— fold vertices are distributed along the external equa- 
tor in the graph of Fig. 10a or form a double ring above 
and below it in the case of the graph Fig. 10b. Toroidal 
deltahedra obtained by embedding the prismatic graphs 
of Fig. 10 on a circular torus are shown in Fig. 11 for 
the case of a 5— fold symmetric toroid with T^ = 35 and 
a 7— fold symmetric toroid with V = 49. In the rest of 
the paper we will reserve the symbol TPn for Dunalp's 
toroids and refer with TPna and TPnb to the other two 
classes of toroids with symmetry group Dnu and unit cell 
of as shown in Fig. 10a and Fig. 10b respectively. 

All defective triangulations presented so far are charac- 
terized by an even number of disclination pairs. Regular 
tessellations of the torus comprising an odd number of 
defects pairs are also possible. Such tessellations are ob- 
tained by combining segments of prismatic and antipris- 
matic lattices with a consequent loss of dihedral symme- 
try. Fig. 12 shows the Voronoi diagram of a toroidal 
lattice, with r — 10/3 and V — 200 vertices, containing 
11 disclination pairs. For an angular length of approxi- 
mately A0 = 7/5 TT the lattice is a prismatic D^d toroid 
while in the remaining 3/5 tt the local structure is that 
of an antiprismatic toroid. The global structure has only 
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FIG. 12: (Color online) Voronoi diagram of a toroidal lattice 
with r = 10/3 and V — 200 vertices. The lattice exhibits 11 
disclination pairs and has Cs symmetry group. 



bilateral symmetry about a sagittal plane dividing the 
lattice in two mirror halves and thus point group Cg- 

In the past few years, alternative constructions of tri- 
angulated tori have been proposed as well as novel ge- 
ometrical and graph-theoretical methods to express the 
coordinates of their three-dimensional structures (see for 
example Kirby^^, Laszlo at aP'^, Diudea et aP^). Here 
we choose to focus on the defect structure associated with 
the two most important class TPn and TAn with groups 
Dnh and Dnd- 



IV. DISCLINATIONS AND SCARS 

A. Isolated Defects Regime 

In this section we analyze the crystalline structures 
arising from the solution of the elastic problem and we 
show how the interplay between the geometry of the em- 
bedding torus, the topology of the lattice, and the me- 
chanical properties of the microscopic units, here encoded 
in the Young modulus Y and the core energy, lead to a 
rich variety of structures whose phase-diagram is pre- 
sented at the end of the section. Since the free energy of 
Eq. (4) is minimized when disclinations best approximate 



the continuum Gaussian curvature of the torus, it is clear 
that disclinations are most likely to be found in regions 
of like-sign Gaussian curvature. Maximal curvature oc- 
curs along the external (positive) and internal (negative) 
equators, which thus constitute preferred regions for the 
appearance of disclinations in the ground state. 

To analyze the elastic free energy (4) we start by con- 
sidering the energies of two opposite sign disclinations 
constrained to lie on the same meridian. The elastic free 
energy of this system is shown in Fig. 14 as a function 
of the angular separation between the two disclinations. 
The energy is minimized for the positive (5— fold) discli- 
nation on the external equator (maximally positive Gaus- 
sian curvature) and the negative (7— fold) disclination on 
the internal equator (maximally negative Gaussian cur- 
vature). The picture emerging from this simple test case 
suggests that a good ansatz for an optimal defect pattern 
is a certain number p of equally spaced -1-1 disclinations 
on the external equator matched by the same number of 
equally spaced —1 disclinations on the internal equator. 
We name this configuration with the symbol Tp, where p 
stands for the total number of disclination pairs. 



T„ 



27rfc 
P 



l<k<p 



TT, ■ 



27rfc 
P 



(26) 



i<fe<p 



where the two pairs of numbers specify the (-0, (f>) co- 
ordinates of the positive and negative disclinations re- 
spectively. A comparison of the energy of different Tp 
configurations, as a function of aspect ratio and disclina- 
tion core energy, is summarized in the phase diagram of 
Fig. 15. We stress here that only Tp configurations with 
p even have an embedding on the torus corresponding 
to lattices of the TP| class. Nevertheless a comparison 
with p— odd configurations can provide additional infor- 
mation on the stability of p— even lattices. For small 
core energies, moreover, thermally excited configurations 
with a large number of defects and similar p— polar dis- 
tributions of topological charge are expected to exhibit 
an elastic energy comparable in magnitude with that of 
these minimal constructions. The defect core energy has 
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FIG. 13: (Color onhne) The screening function Fsitp) for dif- 
ferent values of the aspect ratio r. 
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FIG. 14: (Color online) Elastic energy of a 5 — 7 disclination 
dipole constrained to lie on the same meridian, as a func- 
tion of the angular separation. In the inset, illustration of a 
circular torus of radii Ri > 7?2. Regions of positive and neg- 
ative Gaussian curvature have been shaded in red and blue 
respectively. 



been expressed here in the form: 



2p 



Pc^^c^ qf = 2pec 



(27) 



The core energy ec of a single dischnation depends on the 
details of the crystal-forming material the corresponding 
microscopic interactions. A simple phenomenological ar- 
gument (see for example Ref. 31) gives 



Y 



a 
32^ ' 



where a the lattice spacing. Taking a^ = A/^V, with 
A the area of the torus, yields: 

1 



AY 



10" 



mV^irV 



V 



(28) 



For a system of order V = 10'^ subunits, then, the di- 
mensionless core energy on the left hand side of Eq. (28) 
is of order 10^^. This estimate motivates our choice of 
the scale for ec/iAY) in Fig. 15. 

For dimensionless core energies below 4 • 10^^ and as- 
pect ratios r between 3.68 and 10.12 the ground state 
structure is the TP5 lattice corresponding to a double 
ring of -f 1 and —1 disclinations distributed on the exter- 
nal and internal equators of the torus as the vertices of 
a regular decagon (the Tiq configuration). The TP5 lat- 
tice has dihedral symmetry group D^h- That this struc- 
ture might represent a stable configuration for polygo- 
nal carbon toroids has been conjectured by the authors 
of Ref. 27, based on the argument that the 36° angle 
arising from the insertion of ten pentagonal-heptagonal 
pairs into the lattice would optimize the geometry of 
a nanotorus consistently with the structure of the sp^ 
bonds of the carbon network (unlike the 30° angle of 



the 6— fold symmetric configuration originally proposed 
by Dunlap). In later molecular dynamics simulations, 
Han'^^ found that a 5— fold symmetric lattice, such as 
the one obtained from a (9,0)/(5,5) junction (see Fig. 5), 
is in fact stable for toroids with aspect ratio less then 
r ~ 10. The stability, in this case, results from the strain 
energy per atom being smaller than the binding energy 
of carbon atoms. Irrespective of the direct experimental 
observation of such disclinated toroidal crystals, which is 
still open, we have shown here, from continuum elastic- 
ity, that a 5— fold symmetric lattice indeed constitutes a 
minimum of the elastic energy for a broad range of aspect 
ratios and defect core energies. 

For small aspect ratios the 5— fold symmetric config- 
uration becomes unstable and is replaced by the 9— fold 
symmetric phase Tg. As we mentioned, however, this 
configuration doesn't correspond to a possible triangula- 
tion of the torus. It is likely that the ground sate in this 
regime consists of ten skew disclination pairs as in the an- 
tiprismatic TAn lattice. The latter can be described by 
introducing a further degree of freedom Sip representing 
the angular displacement of defects from the equatorial 
plane: 



TAn: 



i-iy^Si; 



27rfc 



l<k<r 



i-irin-Si.),^-^ 
n 



(29) 



l<fc<r 



A comparison of the TP5 configuration and the TA5 
configuration is shown in Fig. 17 for different values 
of dip. The intersection points of the boundary curves 
with the (5'0— axis has been calculated by extrapolating 



AY 




FIG. 15: (Color online) Phase diagram for Tp configurations 
in the plane {r,e^/AY). For r G [3.68, 10.12] and Ec ~ the 
structure is given by a Tio configuration with symmetry group 
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FIG. 16: (Color online) Two examples of defect free "crop 
circle" toroids with r = 20 and V = 180 (left) and 220 (right). 



the (r, Sijj) data points in the range 5 G [0.07, 0.8] with 
A(JV') = 2.57r • 10-3. For small 5iP and r e [3.3, 7.5] 
the prismatic TP5 configuration is energetically favored. 
For r < 3.3, however, the lattice undergoes a structural 
transition to the TA5 phase. For r > 7.5 the prismatic 
symmetry of the TP5 configuration breaks down again. 
In this regime, however, the elastic energy of both con- 
figurations rapidly rises because of the lower curvature 
and defects disappear. 
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FIG. 17: (Color online) Phase diagram of a 5— fold symmet- 
ric lattice in the plane {r,5ip). For small Sip and r in the 
range [3.3, 7.5] the prismatic TP5 configuration is energeti- 
cally favored. For r < 3.3 the system undergoes a structural 
transition to the antiprismatic phase TA5. 



B. Defect-free Tori 

For aspect ratio r > II the TP5 phase is replaced by 
a defect free configuration (Tq in Fig. 15) so that con- 
figurations with defects are no longer energy minima. 
Any toroidal crystal with aspect ratio larger than ~ 11 is 
than energetically favored to be defect-free. In the thin 
torus limit the ground state structure is directly related 
to the simple problem of finding the most efficient pack- 
ing of congruent equilateral triangles on the torus of a 
given aspect ratio. Given V subunits (vertices) one seeks 
the densest packing of equilateral triangles of edge-length 
a = {A/^ Vy^^ on the torus with aspect ratio r, such 
that each vertex has valence six. Using the planar con- 
struction described in Sec. Ill, the optimal choice of the 
indices (n, rn, I), can be translated into the minimization 
of the following quantity: 



A"'"(r-, V) = n^ + nm + m'^- 



V3 . 



-V, 



(30) 



obtained by equating the magnitude of the chiral vector c 
with that of the sectional circumference of the embedding 
torus, under the constraints: 



I 



V {n+2m:2n+m) 
2 n^+nm+T??.2 



n, m, I Cz 



(31) 



This construction successfully predicts the structure of 
the lattices of Fig. 16. 

So far we have studied the elasticity of toroidal crystals 
exclusively in terms of interacting topological defects on 
a rigid toroidal substrate. Thus the elastic strain due to 
defects and curvature takes the form of pure stretching 
on the tangent plane of the torus and no out-of-plane 
deformation takes place. In a more realistic scenario, a 
crystalline torus would undergo both in-plane stretching 
and out-of-plane bending. The latter implies an energy 
cost: 



F. = ^ 



I (fxH'^ix) =Kb 






(32) 



with H the mean curvature (see Ref. 33). The case of 
defect-free tori is simple enough to incorporate bending 
in the problem and see what the optimal aspect ratio of 
a defect-free torus would be as a function of the Foppl- 
von Karman number 7 = AY/Kf, representing the ratio 
of the stretching energy scale to the bending rigidity. In 
absence of defects the only source of stress is given by 
the curvature. Thus 

Fs^^Jd'xTlix) 

1 + 4r(r2 - 1)5 [1 - log(2 + 2ra)] 



AY 



{- 



2r2 



+ Li2(a2)-2}. (33) 
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FIG. 18: (Color online) Optimal value of the aspect ratio r 
as function of the Foppl-von Karman number 7 = AY/ nb- 
For 7 ~ the Clifford torus with r = \/2 is optimal. Larger 
values of 7 favour instead a "skinnier" torus. 



Summing Eq. (32) and (33) and taking the derivative 
with respect to r (assuming constant area), one obtains 
the foUowing equation for the optimal value of r: 



27r- 



/ 2 



(r2 - 1)5 



l + 2m 



2r(r2 - l)5log(2 + 2ra) 



= . (34) 



The optimal aspect ratio r as obtained form Eq. (34) is 
shown in Fig. 18 as a function of 7. For 7 ~ 0, when 
the major contribution to the elastic energy is given by 
the bending, the optimal geometry is given by the Clif- 
ford torus (r = y/2). If, on the other hand, the in-plane 
stretching dominates, a "skinny" torus (large r) is ener- 
getically favoured. 



C. The Coexistence Regime 

In this section we investigate the possibility of a regime 
of coexistence between isolated disclinations and grain 
boundary "scars" . The existence of scars, first predicted 
in the context of spherical crystallography^^ and later ob- 
served experimentally in spherical droplets coated with 
coUoids^"*'^^ , has become one of the fundamental signa- 
tures of dense geometrically frustrated systems. 

In the regime of large particle numbers, the amount of 
curvature required to screen the stress field of an isolated 
disclination in units of lattice spacing becomes too large 
and disclinations are unstable to grain boundary "scars" 
consisting of a linear array of tightly bound 5 — 7 pairs ra- 
diating from an unpaired disclination^^'^^. In a manifold 
with variable Gaussian curvature this effects leads to a 
regime of coexistence of isolated disclinations (in regions 
of large curvature) and scars. In the case of the torus the 



Gaussian curvature inside (IV'I > '^1'^) is always larger in 
magnitude than that outside (|V'| < 7''/2) for any aspect 
ratio and so we may expect a regime in which the neg- 
ative internal curvature is still large enough to support 
the existence of isolated 7— fold disclinations, while on 
the exterior of the torus disclinations are delocalized in 
the form of positively charged grain boundary scars. 



To check this hypothesis we compare the energy of the 
TP5 lattice previously described with that of "scarred" 
configurations obtained by decorating the original toroid 
in such a way that each -1-1 disclination on the exter- 
nal equator is replaced by a 5 — 7 — 5 mini-scar. The 
result of this comparison is summarized in the phase dia- 
gram of Fig. 19 in terms of r and the number of vertices 
of the triangular lattice V (the corresponding hexagonal 
lattice has twice the number of vertices, i.e. Vhex — '2V). 

V can be derived from the angular separation of neigh- 
boring disclinations in the same scar by approximating 

V ~ A/Av, with Ay = ^c? the area of a hexagonal 
Voronoi cell and a the lattice spacing. When the aspect 
ratio is increased from 1 to 6.8 the range of the curva- 
ture screening becomes shorter and the number of sub- 
units required to destroy the stability of the TP5 lattice 
decreases. For r > 6.8, however, the geodesic distance 
between the two equators of the torus becomes too small 
and the repulsion between like-sign defects takes over. 
Thus the trend is inverted. 
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FIG. 19: (Color online) Isolated defects and scar phases in 
the (r, 1/) plane. When the number of vertices V increases 
the range of the screening curvature becomes smaller than 
one lattice spacing and disclinations appear delocalized in the 
form of a 5 — 7 — 5 grain boundary mini-scar. 
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V. NUMERICAL EXPERIMENTS 

In this section we report the resuh of a numerical mini- 
mization of a system of V point-hke particles constrained 
to he on the surface of a torus and interacting via a pair 
potential of the form Uij ~ l/\xi — Xj\^ where | • | de- 
notes the Euclidean distance in M.^. The problem of find- 
ing the minimal energy configuration of repulsively in- 
teracting points on a 2-manifold has become a standard 
problem of potential theory and has its paradigm in the 
classical Thomson problem on the sphere. The choice 
of the cubic potential is motivated here by the so called 
"poppy seed bagel theorem" ^°, according to which the 
configuration of points that minimizes the Riesz energy 
E ~ X]i<7 1/l^i ^ ^j\^ on a rectifiable manifold of Haus- 
dorff dimension d is uniformly distributed on the mani- 
fold for s > d. In the case of a torus of revolution this 
implies that for small s the points are mostly distributed 
on the exterior of the torus (the interior becomes com- 
pletely empty in the limit V^ ^ oo). As s is increased, 
however, the points cover a progressively larger portion of 
the surface. The distribution becomes uniform for s > 2. 
On the other hand, since the number of local minima of 
the Riesz energy increases with s, it is practical to choose 
a value not much larger than two. The choice s = 3 has 
the further advantage of modelling a real physical system 
of neutral colloidal particles assembled at an interface''^ 
and is therefore suitable for direct comparison with ex- 
periments on colloidal suspensions. 

To construct low energy configurations we adopt an 
carefully designed hybrid optimization algorithm named 
Tapping (TA). Like other hybrid algorithms, TA consists 
in a combination of fast local optimizations and global 
stochastic moves designed to release the system from the 
local minimum to which it is confined at the end of a 
local minimization step. A more detailed description of 
our algorithm is reported in Appendix C. We study four 
different aspect ratios: r = 3,4,6 and 20. For each as- 
pect ratio we consider several different particle numbers 
up to y = 1000 and each simulation is performed for 10^ 
to 10^ TA iterations. 

The lowest energies found, as well as the number of de- 
fects in the corresponding configuration, are reported in 
Table II for a selected set of systems. The corresponding 
lattices are shown in Fig. 20 

The lattices are best presented using a Voronoi con- 
struction corresponding to the dual lattice of the Delau- 
nay triangulation. Here pentagonal faces are colored in 
red while heptagonal faces are colored in blue. The com- 
plete set of data produced in our simulations together 
with a collection of interactive 3D graphics for each low 
energy configuration found is available on-line^^. For 
fewer than y ~ 180 particles the results of our numer- 
ical minimization are in good agreement with the con- 
tinuum elastic theory. In particular for 180 < V < 500 
and r = 3, 4 and 6, we always find minimal energy con- 
figurations consisting of ten 5— fold disclinations on the 
outside of the torus and ten 7— fold disclinations in the 



inside as predicted by the elastic theory in the regime of 
ec/{AY) ~ 0. For r = 20 and F > 110 we also find the 
lowest energy configurations to be defect free. 

For small numbers of particles we don't expect the 
continuum approximation to accurately describe the low- 
est energy structure of the toroidal clusters presented in 
Fig. 20. Loosely speaking the limit of validity of the 
elastic theory can be quantified by requiring the aver- 
age lattice spacing a = 2'k[RiR2/ {^V)]^^^ to be much 
smaller than the radius i?2 of the torus. This condition 
requires V to be of order 500 particles for a torus with 
aspect ratio r = 3. Remarkably, good agreement between 
the theory and simulations is found starting from much 
smaller values of V and in some cases (see the following 
discussion on the configuration with r = 3 and V = 130), 
we already observe the onset of the ideal behavior pre- 
dicted by theory for a ^ i?2 . The occurrence of a ground 
state configuration with exact prismatic or antiprismatic 
symmetry, in particular, is only possible when the num- 
ber of particles V belongs to a specific sequence of "magic 
numbers" described in Sec. III. Nevertheless for V out- 
side such a sequence it is still possible to observe in the 
ground state a predominant prismatic or antiprismatic 
character depending on the aspect ratio. 



V 


rn,in ± 0.05 


rmax ± 0.05 


16 


1.0 


1.6 


20 


1.4 


2.6 


24 


1.8 


3.4 


28 


2.4 


4.0 


32 


2.9 


4.6 



TABLE I: Maximum and minimum aspect ratio for which the 
toroidal antiprisms are a global minimum. 

Some configurations deserve special attention. For 

V = 16, 20, 24, 28 and 32 and r within a specific range 
(see Table I) the global minima are represented by the 
second to sixth toroidal antiprims discussed in Sec. III. 
The drilled icosahedron, on the other hand, would re- 
quire the aspect ratio be less than one, as can be under- 
stood from Table I, and is therefore never a minimum for 

V = 12. We next describe the salient features of the four 
aspect ratios simulated. 

r = 3) The smallest minimal energy state with D^h 
symmetry is obtained for V ~ 35. It features ten discli- 
nation pairs and belongs to the class of TP5a graphs. 
For V — A2, a. TP6b lattice is obtained with no defects 
along the two equators. Two 6— fold chiral configurations 
are obtained for T^ = 60 and 126. The global minimum 
obtained for V — ViQ displays a fascinating example of 
5— fold antiprismatic symmetry with the ten isolated neg- 
ative disclinations in the interior of the torus replaced by 
a simplicial complex consisting of five triangles with a 
common 5— fold apex and four 7— fold coordinated ver- 
tices along the base. A peculiar example is also rep- 
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resented by the minimum obtained for V = 180. The 
lattice exhibits the typical pattern of a TPn graph with 
(3, 1,4, 7)— type unit cell. The angular distance between 
neighboring disclinations is S(j) ~ 27r/9. Since a pris- 
matic graph cannot have an odd number of disclination 
pairs the toroidal lattice is closed by a simplicial complex 
consisting of two positive disclinations on the exterior of 
the torus at the opposite sides of the external equator 
and two negative disclinations in the interior arranged 
similarly. The total number of disclination pairs is there- 
fore ten. The typical pattern of the 5— fold antiprismatic 
toroid can be found in all configurations with V > 200. 
A single 5 — 7—5 scar appears in the r — 3 configura- 
tions at V^ = 420, while larger lattices (i.e V = 460 and 
500) also feature 4— fold disclinations in the interior of 
the torus. It is not clear, however, whether the presence 
of disclinations with topological charge |g| > 1 is a gen- 
uine property of the ground state or rather an artifact 
due to a misconvergencc of our algorithm. 

r = 4) An interesting feature is observed at t^ = 42. 
As in the case of r = 3 we also find a minimum with 
Dg/i symmetry group, but unlike the latter configura- 
tion, it belongs to the TP6a class and has 5— fold discli- 
nations along the external equator. A TP7b configura- 
tion is obtained again for V — 49. For F = 66 and 

V = 104 the global minimum is achieved by two spec- 
tacular antiprismatic configurations with Dud and I?i3d 
symmetry group respectively. These toroids can be ob- 
tained from the toroidal antiprisms discussed in Sec. Ill 
by splitting^^ one or more times the initial set of 5— fold 
vertices. Thus starting from a 11— fold toroidal antiprism 
with l/ = 4xll = 44 vertices and spfitting all 22 5-fold 
vertices one obtains y = 44 -I- 22 = 66 vertices. Split- 
ting twice all 26 5— fold vertices of a 13— fold toroidal 
antiprism with ]/ = 4 x 13 = 52, on the other hand, we 
have F = 52 + 2 X 26 = 104. For V = 120 the global 
minimum is represented by a fascinating lattice of TP5a 
type. Lattices with V = 121 , 125 and 126 resemble very 
closely the structure of a TP5 graph while for V ~ 260 
the lattice has a more antiprismatic character with ten 
defect pairs. 

r — 6) Three defect-free configurations are found at 

V = 87, 112 and 116. The case F = 87 is a particu- 
lar example of a defect-free lattice that cannot be ob- 
tained form the planar construction reviewed in Sec. III. 
It consists of 29 octahedra connected in the form of a 
chain. Since each octahedron is attached to other two, 
it contributes with six faces to the total face count. 
Thus F = 6 X 29 == 174 and y = 174/2 = 87. For 
180 < V < 460 we always find configurations with ten 
disclination pairs as expected from continuum elasticity. 
For V > 460 the regime of coexistence between isolated 
disclinations and scars described in Sec. IV C is ob- 
served. The delocalization of isolated disclinations into 
scars, however, doesn't take place at each defective site si- 
multaneously and the regime of coexistence between pos- 
itively charged scars and isolated 7— fold disclinations is 
preceded by a phase with isolated 5— and 7— fold discli- 
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20 
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222 
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102100.926892 
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1 
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12 





146139.605664 
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10 
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10 





199812.441922 
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11 


398 


11 





339147.966681 
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14 


426 


18 





425754.968401 
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2 


16 


462 


20 





524508.172150 
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17 


963 


19 





2965940.674307 
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22 


22 


22 





4905.964854 
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26 


52 
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15598.534409 
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15 


75 


15 





15984.990289 




113 





19 


75 


19 





19237.981548 




117 





22 


73 


22 





21007.172188 
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12 


95 


12 





21914.283713 
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120 





10 


100 


10 





22371.402771 




121 





12 


97 


12 





22859.735385 




125 





10 


105 


10 





24816.591295 




126 





10 


106 


10 





25311.298095 




180 





10 


160 


10 





62142.129092 




220 





10 


200 


10 





102919.127703 




260 





10 


240 


10 





156499.285669 




300 





10 


280 


10 





223997.341297 




340 





10 


320 


10 





306568.539431 




420 





13 


394 


13 





520431.653442 




460 





11 


438 


11 





653485.181907 




500 





14 


472 


14 





805206.972227 




87 








87 








17765.124942 




108 





12 


84 


12 





30894.374674 




112 








112 








33902.717714 
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115 








36254.709031 


6 


116 








116 








37074.949162 




180 





10 


160 


10 





112810.451302 




220 





10 


200 


10 





187146.462245 




260 





10 


240 


10 





284907.016076 




340 





10 


320 


10 





559161.546358 




420 





10 


400 


10 





950488.931696 




500 





13 


474 


13 





1471923.063515 




1000 





30 


940 


30 





8351619.696538 




160 








160 








463967.242489 




170 








160 








543799.839326 


20 


180 








180 








631751.371902 




220 








220 








1065625.748639 




260 








260 








1636942.532923 




300 








300 








2370110.403872 



TABLE II: Low energy configuration for a selected number 
of toroidal lattices with aspect ratios r = 3, 4, 6 and 20. For 
each aspect ratio the table displays the number of particles 
V, the lowest energy found and the number of fc— fold vertices 
T4 with fe = 4-8. 
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nations and scars. For tori with aspect ratio as large as 
r = 20 we find defect-free ground states every time it 
is possible to construct a purely 6— valent toroidal graph 
with the same number of vertices V . 



VI. THE FAT TORUS LIMIT 

We have seen that disclination defects, forbidden in 
the lowest energy state of a planar crystal, may be ener- 
getically favored on a substrate of non-vanishing Gaus- 
sian curvature. It is therefore natural to ask whether 
large curvature can completely destroy crystalline or- 
der by driving the proliferation of a sufhciently high 
density of defects. The resulting state would be amor- 
phous. The problem of generating amorphous structures 
by tiling a two-dimensional curved space with identical 
rigid subunits has drawn attention over the years, par- 
ticularly through the connection to the structure of such 
disordered materials as supercooled liquids and metallic 
glasses. Since the work of Frank'^'' the notion of geometri- 
cal frustration arises frequently in investigations of super- 
cooled liquids and the glass transition. A paradigmatic 
example is represented by the icosahedral order in metal- 
lic liquids and glasses which, although locally favored, 
cannot propagate throughout all of three-dimensional 
Euclidean space. A two-dimensional analog, consisting 
of a liquid of monodisperse hard disks in a 2-manifold 
of constant negative Gaussian curvature (the hyperbolic 
plane) was first proposed by Nelson and coworkers in 
1983'^*. In such a system the impossibility of covering 
the entire manifold with a 6-fold coordinated array of 
disks mimics many aspects of the geometrical frustration 
of icosahedral order in three dimensions. In all these 
models of geometrical frustration, however, the origin of 
the disorder is primarily due to the short-range nature of 
the potential between the subunits. In a more realistic 
setting, part of the frustration is relieved by the fact that 
hexagonal unit cells can compress in order to match the 
underlying geometry. 

The embedding of a triangular lattice on an axisym- 
metric torus, provides a particularly suitable playground 
to study curvature-driven disorder. When r ^ 1 the 
Gaussian curvature on the inside of the torus grows like 
l/(r — 1) and diverges on the internal equator &t ip = i:. 
We thus expect a high density of defects in the vicinity 
of the curvature singularity and a resultant loss of the 
local 6— fold bond orientational order. In this regime the 
system will have crystalline regions on the outside of the 
torus and amorphous regions near the curvature singu- 
larity. 

In this section we substantiate this claim analytically 
based on the elastic theory of continuous distributions 
of edge dislocations on a "fat" torus. Our argument is 
based on the following construction. As a consequence 
of the curvature singularity the surface area of an arbi- 
trary wedge of angular width /S.(j) becomes smaller and 
smaller as the sectional angle f/' increases and vanishes 



at V' = TT. If a defect-free lattice is embedded on such 
a wedge, Bragg rows will become closer and closer as 
the singularity is approached with a consequent rise in 
the elastic energy (see Fig. 21). An intuitive way to re- 
duce the distortion of the lattice is to recursively remove 
Bragg rows as one approaches the point if) = n (see Fig. 
22). This is equivalent to introducing a growing den- 
sity of edge dislocations. This dislocation "cloud" will 
ultimately disorder the system by destroying the local 
6— fold bond orientational order. One might therefore 
view the curvature as playing the role of a local effec- 
tive temperature which can drive "melting" by liberating 
disclinations and dislocations. In two-dimensional non- 
Euclidean crystals at T = 0, however, the mechanism for 
dislocation proliferation is fundamentally different from 
the usual thermal melting. While the latter is governed 
by an entropy gain due to unbinding of dislocation pairs, 
the amorphization at T = is due to the adjustment 
of the lattice to the geometry of the embedding mani- 
fold via the proliferation of defects and the consequent 
release of elastic stress. A similar phenomenon occurs 
in the disorder-driven amorphization of vortex lattices in 
type-II and high- Tc superconductors^®. 

Since the shrinking area per plaquette on the inside 
of the torus necessitates a high density of dislocations 
we may approximate the dislocation cloud in this region 
by a continuous distribution of Burgers vector density 
b. Minimizing the elastic energy with respect to b yields 
a variational equation from which the optimal disloca- 
tion density can be calculated as a function of the ra- 
tio ed/iYR^) between the dislocation core energy e^ and 
elastic energy scale YB? with R= Ri = i?2 • 

As a starting point, we calculate the Green function 
GL{x,y) in the fat torus limit r ^> 1 (i.e. k ^ oo and 
u; — > 0, see Eq. (13)). The conformal angle ^ in this limit 



is: 



lim £ 



tan 



V' 



and to leading order of k. we have: 



167r2 



i'-'i 



IGtt 



„ lb tan — , 

2 ^ 47r 2 ' 



-^Re{Li2(ae*'^)} 

47r^ 



I^^ +^log(cos- 



To handle the limit of the Jacobi theta function we can 



take u = IS.zJK^ q 



— e » and calculate the limit 



(J ^ 1. This can be done by using the modular transfor- 
mation properties of Jacobi functions'^" : 



^i(7|-7)=-«(-*T)5e^7?i(w|r) 



(35) 



Thus t' = -1/r = iK/2, u' = u/t = Az/(2i) and q' = 
e" — e^^' , where: 



lim -di (u, q) = lim i \ — e^^^'di (u ,q) 
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FIG. 20: (Color online) Selected low energy configurations for toroidal lattices of aspect ratio r = 3, 4, 6 and 20. Lattices are 
labeled by (r, V), with V the number of particles. 
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FIG. 21: (Color online) Top view of a defect free triangulation 
of a fat torus with [n,m,l) = (10, 10,20) and V = 400. The 
corresponding elastic energy becomes very high in the interior 
of the torus where the triangles are more compressed to match 
the reduction of surface area. 




FIG. 22: (Color online) A schematic example of a dislocation 
pile-up on a square lattice resulting from the shrinking of the 
area on a regular wedge of a fat torus. 



This is easily evaluated by means of the expansion: 



■di {u,q) — 2(7 4 sin u + o Iq-^j 



Taking the logarithm and neglecting irrelevant constant 
terms, we obtain: 



log 



^1 



log 



sinh 



which finally leads to: 



GL{lp,(l),tp',<l>') 



--—J tan — 
47r2 2 



1 

2^ 



log 



sinh 



(36) 



with z ~ tan('i/;/2)+i0. With the Green function in hand, 
we can calculate the effect of the curvature singularity at 
-0 = TT on the distribution of defects. Let b be the Burgers 
vector density of the dislocation cloud. Hereafter we work 
in a local frame, so that: 



b = b^g,, + b^g^ , 



(37) 



with gi = diR a, basis vector in the tangent plane of the 
torus whose points are specified by the three-dimensional 
Euclidean vector R. The quantity b has to be such that: 



(fix b{x) — bo 



with bo the total Burger's vector in a generic domain 
D. Because on a closed manifold dislocation lines cannot 
terminate on the boundary, extending the integration to 
the whole torus we have: 



drxb{x) — . 



(38) 



Since the basis vectors gi in Eq. (37) have the dimension 
of length, contravariant coordinates V have dimensions 
of an inverse area. Assuming all defects to be paired in 
the form of dislocations (i.e. qi = Q everywhere), the 
total energy of the crystal reads: 



F = ^ J fixT^x) + Cd J fix\b{x)\' 
where e^ is the dislocation core energy and 



(39) 



Ibl- 



gijb'fy' = g^^{b'^f +g^^ 



The function r(a;) encoding the elastic stress due to the 
curvature and the screening contribution of the disloca- 
tion cloud obeys 



-Agr{x) = elS/,b'^{x)-K{x), 



(40) 



where e^ is the Levi-Civita antisymmetric tensor on the 
torus: 



eV0 



-£0-0 — \/g , 



4^9^ke'' 



The stress function r(a;) can be expressed in the form 

r(a;)-rrf(^,0)-r,(0) with 



Y 



log 



2(1 + cost/") 



Y 



fiyelW,b'^{y)GL{x,y). 



(41a) 
(41b) 
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Taking advantage of the closeness of the torus we can 
integrate Eq. (41b) by parts so that: 



Y 



= -Jd'yelb''iy)d,GL{x,y). (42) 



Now we want reduce the integral term in Eq. (42) to a 
more friendly functional of b, suitable for a variational 
approach. Given the azimuthal symmetry we assume 
that all dislocations are aligned along b = h'^g^,. Even 
though not necessarily true, we argue this to be a rea- 
sonable work hypothesis as well as a solid starting point 
to capture the essential physics of the fat limit. In this 
case Tdi'^l^, <j}) = ^di^-') can be recast in the form 



Y 2tt 



[ip' + sin tJj' + n sgn{ip - V')] ■ (43) 



Substituting Eq. (43) and (41a) in Eq. (39) and minimiz- 
ing with respect to b'^ we can now write the variational 
equation: 



4edi?2(l + cosV')'&'^(^)+ / d^^Td{4'')sgni^-,p') 
d^V5r.(^')sgn(V'-V'')- (44) 



By inverting the order of integration in the integral on 
the right hand side, Eq (44) can be expressed in the form 
of a Fredholm equation of the second kind: 

XB{^)- f d^'B{ij')lC{i,,i,')^f{^), (45) 



where A = ed/lYR^), B{i;) ^ R^{1 + cosV')^6'^(^) and 
the kernel IC{ijj,ip') is given by: 



+ |^_^'| + 2cosi±i:^sinM_il|. (46) 



The function f{tp) on the right hand side of Eq. (45) is 
given by: 

fi^) = -\j d^'(i + cosV')r.(V') 

= i{[log2(l + cos^) - 2] siuT/; - 2Cl2(V^ + tt)} (47) 

where CI2 is the Clausen function (see Ref. 43, pp. 1005- 
1006) defined as: 



Cl2(x) 



/ dx log I 2 sin - I = > 
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FIG. 23: (Color online) The Burgers vector component &(■(/)) 
for different choices of A. 



As previously noted the dislocation core energy is ed is 
much smaller than the elastic energy scale YB?'. Eq (45) 
is then suitable to be solved in powers of the dimension- 
less number A: 

B{ij) - Bo(V') + \Biii,) + X^B^i^) + • • • 

The corrections to the zero-order term Bo{ip) can be cal- 
culated recursively by solving a set of Fredholm equations 
of the first kind: 

Bk-i{i^)= f rf^'Bfe(7A)/C(^,V') k>l. 

The function Bo{ip) associated with the Burgers vector 
density of the dislocation cloud in the limit A — > 0, on 
the other hand, can be calculated directly from Eq. (40) 
by setting the effective topological charge density on the 
right hand side to zero: 



elV^b''{x)-K{x)=0. 



(48) 



For a torus of revolution the only nonzero Christoffel 
symbols are 






R2 sin ijj 



Ri -f i?2 cos ^ ' 

r^0 - R2^ sin ^(i?! + i?2 COS V) . 

Since fe''^ = by assumption, the first term in Eq. (48) 
can be expressed as: 



^^V.b^ = el{d4 



■ ^^i 






= (1 -h COS V')5v,6'^ - 2 sin V' fe"^ , 
and Eq. (48) becomes an ordinary differential equation 

2 sin tp , cos V' 



d^ 



I + cost/i R^IX + cosijjY 



(49) 
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whose solution is given by 



sin^ 



i?2(l + cosV')2 



(50) 



so that Bo{ip) = sin-0. The Burgers vector density b'^ 
obtained by a numerical solution of Eq. (45) is shown 
in Fig. 23 for different values of A. The Burgers vector 
density is measured in units of R~^. The function b'^ has 
cubic singularities at ip — ±7r and is approximately zero 
on the outside of the torus. The solid blue curve in Fig. 23 
represents the zeroth order solution of Eq.(50). 

Now, in the theory of dislocation mediated melting a 
system at the solid liquid phase boundary is described as 
a crystalline solid saturated with dislocations. In three- 
dimensions, in particular, there is a strong experimental 
evidence of the existence of a critical dislocation density 
at the melting point p{Tm) ~ 0.66"^ where b is the length 
of the length of the smallest perfect-dislocation Burgers 
vector*''. Several theoretical works have motivated this 
evidence both for three-dimensional solids and vortex lat- 
tices in super conductors'^^. On the other hand, given the 
existence of such a critical density, its value can be em- 
pirically used to determine whether a system is in a solid 
or liquid-like phase in the same spirit as the Lindemann 
criterion. With this goal in mind we can calculate the 
dislocation density by requiring |6| ^ pa with p the den- 
sity of single lattice spacing dislocations. This yields: 



p(^, V)a^ = 2tt 



'V3 



V 



tan 



i, 



+ o(A) (51) 



Solving p[xl}, V)a'^ — 0.6 as a function of i/' and V we ob- 
tain the diagram of Fig. 24. As expected the inside of the 
torus contains an amorphous region whose angular size 
decreases with the number of vertices T^ as a consequence 
of the reduction of the lattice spacing. 
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FIG. 24: (Color online) Phase diagram for curvature driven 
amorphization. The inside of the torus contains an amor- 
phous region whose angular size decreases with the number 
of vertices V as a. consequence of the reduction of the lattice 
spacing. 



7— fold disclinations in the inside of the torus together 
with scars on the outside. For a torus of aspect ratio 
close to one we showed how a diverging Gaussian cur- 
vature on the internal equator is responsible for the re- 
markable occurrence of a curvature-driven transition of 
the system to a disordered, liquid-like, state. The predic- 
tions of our elastic theory were compared with the results 
of a numerical study of a system of point-like particles 
constrained on the surface of a torus and interacting via 
a short range potential, with good agreement. From a 
purely geometrical point of view we have introduced a 
number of novel toroidal polyhedra as well as a new con- 
struction and classification scheme for certain types of 
prismatic tori. 



VII. DISCUSSION AND CONCLUSIONS 



In this article we reported a comprehensive analysis 
of the ground state properties of torodial crystals. Us- 
ing the elastic theory of defects on curved substrates we 
identified the ground state structures of an arbitrary crys- 
talline torus as a function of the aspect ratio r and the 
ratio €c/{AY) of the defect core energy to the elastic 
energy scale. We showed that for a large range of as- 
pect ratios and core energies the minimal energy struc- 
ture of a toroidal crystal has ten disclinations pairs and 
symmetry group -D5/1, as first conjectured by Lambin et 
al more than ten years ago^^. For large system sizes 
we proved isolated disclinations are unstable to grain- 
boundary scars consisting of chains of tightly bound 5 — 7 
pairs radiating form an unpaired disclination. On a torus, 
where the Gaussian curvature on the inside is always 
larger in magnitude than that on the outside, the oc- 
currence of scars is marked by a state featuring isolated 
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APPENDIX A: THE GREEN FUNCTION ON 
THE TORUS 



where A„ and fXm are given by: 



The modified Laplacian Green function on a periodic 
rectangle of edges pi and p2 can be conveniently calcu- 
lated in the form: 






u\{x)u),{y) 
X 



(Al) 



where u\ is the eigenfunction of the Laplace operator 
with periodic boundary conditions: 



27rn 2TTm 

A„ = Mm = n, TO = 0, ±1, ±2. 

Pi P2 



and the eigenvalue A is given by: 



^ — ^K ^ /"ri 



(A4) 



Aux{x) = \ux{x). 



(A2) 



such that: 



Calling for simplicity x — (x, y) and y — [S,, ?y), the func- 
tion Go is given by: 



ua(C,0) =ux{£.,P2) 

In Cartesian coordinates the eigenfunctions arc simple 
plane waves of the form: 



u\{^,v) = 



= «(A„C + Mm»?) 



VpTp2 



(A3) 



Go{x,y) 



PlP2 



(n,m)^(0,0) 






(A5) 



Although (A5) is simple, it is very useful to rewrite it in 
terms of elliptic functions. Noting that the odd terms in 
(A5) cancel we have: 



Go{x,y) = -- 



PlP2 



E 



(n,m)#(0,0) 



COS A„ (x ~ ^) COS M„i jy-i]) 

^n +/^m 



PlP2 



°2, cos^^^{y-'n) 



m—l 



f 27Trn\' 
\ P2 J 



OO 

E 

n— 1 m— — oo 



E 



cos ^(x - COS ^{y-ri) 

( 2-K71 \ \ ( 2TT'm \ 

PI ) ^ \ P2 ) 



(A6) 



r 



An equivalent expression can be obtained by isolating The second sum in Eq. (A6) can be evaluated with the 
the jTi = contribution in the sum rather than the n = help of the Poisson summation formula: 
term. The first sum in Eq. (A6) can be evaluated easily 
by using: 



E 



COSfcx TT Trlxl X 



k=l 

Thus we have: 



/c2 
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(A7) 



2 °2, cos22™i(j/-ry) 
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2 

-^-P2\y-'n\ + \y-v? 



(A8) 



OO <.oo 



^J^ -^ /'OO 

\J f{m)cospm— 2_. / dt f (t) cos{2k'K + p) t . 

(A9) 



r?i— — OO 



fc=-c 



In particular, if we choose: 



p=—{y-v). fim) = - — -2 — ^ 

P2 I 27m \ I / 27rm\ 

V PI / \ P2 J 



we can write the second sum in Eq. (A6) as: 



K{x~^,y-T]) 



2 ^ 2^n, ,, ,^ cos^iy-fj) 

> cos (X — £,) > cos 
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(AlO) 



The integral can be easily calculated by considering: 



dx 



cos LUX TT 



a^ + x'^ 2a 



Thus Eq. (AlO) becomes: 
K{x- ^,y-r]) - 



1 ^ ^ ^~^\P2k+y-v\ 

2^ ^ ^ 



k— — oo n— 1 



cos— (x-e) (All) 
Pi 



The sum in n can be calculated by noting: 
cos{2TTny) — — log ll — e^'^l 

r) 



(A12) 



where a = 27r(a; ± iy) with x > and arbitrary choice 
of sign. Thus separating positive and negative values of 
k in Eq. (All) and using Eq. (A12) the function K{x — 
^,y — 1]) takes the final form: 



2TrK{x -^,y-v) ^ log 



+ Elog 



fe=i 



1 — e PI 



27r, 



(2-0 



l-2g^''cos— (z-o+r'' 
Pi 



(A13) 



where: 



z = x + iy 
C^^ + iV 



9 = e 



The second term in Eq. (A13) can be expressed in terms 
of the Jacobi theta function {)i{u,q) defined as: 

00 
^1 (u, q) = 2qismuY[{l- 2g2" cos 2u + q^") (l - q^") . 

n=l 

Another useful relation can be obtain by taking the 
derivative of i?i('u, q) with respect to u: 

l^ tltlA ^ lin, ^itlA ^ ^;(o, ,) . (A14) 

u^o sin u u^Q cos u 

Thus we have: 
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2gi^(l-g2fe)^ (A15) 



In this way we can express: 

00 
Y[{1^2q^''cos2u + q*'') = 



(2(74)3 siuM 



fe=i 
Then taking: 

Pi 
and substituting in Eq. (A13) we obtain: 
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(A17) 



Combining Eq. (A17) with Eq. (A8) we conclude that: 

Goix, y) = H{y - ry) + K{x -^,y-rj) 



log2 1 , ,2 1 , 
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An alternative notation frequently used for the Jacobi 
theta function is: 

i}i{u\T)^d{u,q) q^e'""^. 

With this choice we can write the Green function in the 
final form: 
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k=l 



APPENDIX B: DERIVATION OF THE 
FUNCTIONS r,(a;) AND rd{x,Xk) 

In this appendix we derive the analytical expression 
for the stress functions Ts{x) and Td{x,Xk) in Eq. (18). 
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The former is given by the integral: 

r.(a^) 



Y 



= J d^yK{y)[Go{x,y) - {Go{;y))] (Bl) 



It is convenient to keep the Green function Go{x,y) in 
the form of Eq. (A6). Thus we have: 



h^ Jd^yKiy)Go{x,y) 



PI r ,,, „^ cos^(g-0 
— / dV'cos^^ -, 



n=l 



which using Eq. (A7) and carrying out the integrals 
yields: 



h = log 



2(r2 - 1) 



(r + cos tp){r + \/r^ — 1) 



(B2) 



To calculate the second integral in Eq. (Bl) it is conve- 
nient to invert the order of integration and use the result 
of Eq. (B2): 

I, = I d'yK{y){G,{;y)) = I ^TsAy) 
4(r2 - 1) 



= log 



-V^^2— J 



Combining Eq. (B2) and (B3) we obtain: 

Ts{x 
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= log 
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2(r + cos V') 
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The elastic stress produced at the point a; by a discli- 
nation at Xk is given by the function: 



rd(x,Xk) 

Y 



Go{x,Xk) - {Go{-,Xk)) 



(B5) 



To calculated the average {Go{-,x)) one can again start 
from Go{x,y) in the form of a series and use Eq. (A7). 
This yields 



A{Go{-,x))=-^R^R2 + l 
2 



-^J\^^'^9\i-i'?- (B6) 



To calculate the integrals in Eq. (B6) it is convenient to 
expand the conformal angle ^ in Fourier harmonics: 



with: 



C(V') = 2J ^" ^^'^ ^'^ 



K - (-1)"] 



(B7) 



where a = \/r'^ — 1 — r. Thus we have: 
h= r dyj'V9\C-e\ 

J— TT 



= 2i?ii?2^e + 2i?^log 
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with 



dV'' ^ — yj — (cos n^ — cos rnr) 



-2 „;,2 



= -K(7r^ - ?A^) - KLi2(-a) + KRc{Li2(ae'''')} (B9) 

where Re{-} stands for the real part and 

00 
Li2(^) = E 



00 9 



is the usual Eulcr's dilogarithm (see Ref. 43 pp 1004- 
1005). The second integral in Eq. (B6) is give by: 

= RiR2 di^'\£,-^'\^-2TTKRl\og ^ 2_L . 
J —TT \_T ^ y T 1 

(BIO) 

To calculate the integral in Eq. (B) one uses Parseval's 
identity: 

where f{x) an arbitrary square-integrable function on the 
interval [— tt, tt] with Fourier series: 

00 
f{x) = — + 2^{0'n cosnx + bn sinnx) . 

n=l 

Thus we have: 
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Using Eq. (B), (B) and (Bll) we conclude that: 
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APPENDIX C: CLUSTER OPTIMIZATION VIA 
TAPPING 

The Tapping Algorithm (TA) is a hybrid algorithm 
designed to find the optimal crystalline structure of par- 
ticle systems constrained to lie on a curved surface and 
interacting with a long range potential of the form Uij = 
l/|ri — Tjl*. Hybrid algorithms, such as Basin-Hopping''''' 
and Minima-Hopping^^ , have been successfully employed 
throughout the years to predict the crystalline structure 
of molecular clusters and proteins. In general they com- 
bine fast local minimizations with global moves whose 
goal is to release the system from the local minimum it 
is confined at the end of a local minimization step. 

A typical hybrid optimization routine can be summa- 
rized in the following two steps: 1) after all the inde- 
pendent variables have been randomly initialized, a local 
optimization is performed and a local minimum x is de- 
termined; 2) from X a new configuration y is constructed 
by applying a global (generally stochastic) move. The 
new configuration y is then used as starting point for a 
new local optimization step. The two steps are iterated 
until a stopping criterion is satisfied. The goal is thus to 
explore the largest possible number of local minima and 
avoid visiting the same minimum too often. 

The crucial point in designing an effective hybrid algo- 
rithm is clearly the choice of the global move. There is 
no general rule to identify a successful global transforma- 
tion X ^ y and physical intuition and prior experience 
are typically the only guidelines. In the case of Basin- 
Hopping, for instance, the global transformation consists 
in a Monte Carlo move in which all the particles of the 
system are randomly displaced in order to construct a 
new initial configuration from which a new trial mini- 
mum is obtained. The step is accepted with probability 
exp(— /3Ay), where AV is the energy difference between 
the new and previous minimum and P is an inverse tem- 
perature adjusted to obtain a 50% acceptance ratio. In 
the case of Minima- Hopping the escape step is performed 
by a short Molecular Dynamics simulation by assigning 
the particles a fixed kinetic energy. 

The global move adopted in TA is inspired by the pro- 
cess of close packing of spherical objects by tapping and 
is motivated by the well established role of topological 
defects in determining the order of two-dimensional non- 
Euclidean crystals as well as the picture of the potential 
energy surface (PES) of such systems as a multi-funnel 
landscape. Consider a system of say spherical objects 
confined in a two-dimensional box with an initial disor- 
dered configuration. A common way to bias the system 
toward a close-packed configuration is to provide it ki- 
netic energy by gently tapping the box. If the system is 
populated by locally ordered regions (i.e. grains) sepa- 
rated by clusters of defects, the primary effect of tapping 
is to produce a glide of defects inside the crystals with 
a subsequent rearrangement of grains. This mechanism 
can be reproduced numerically in the following way. The 
algorithm starts with a random distribution of particles 



and rapidly quenches the system by performing a fast 
local minimization. Once particles are trapped in a lo- 
cal minimum, defects are identified by a Delaunay tri- 
angulation of the lattice. Then the system is tapped by 
adding to the defect positions a random displacement. 
The magnitude of the displacement is given by the typ- 
ical spacing associated with the particles number times 
a factor A which represents the tapping strength. This 
factor is initially set to 10^'^. After defects have been 
moved a new local minimization is performed in order to 
construct the trial configuration y. The energy of this 
configuration is compared with the energy of the previ- 
ous minimum and the move is accepted if their difference 
is larger than some tolerance factor ee- If, on the other 
hand, the energy difference is smaller than e^, the sys- 
tem has relaxed again to the same minimum. In this case 
the tapping strength is increased of a factor 10 and the 
process is repeated until the system successfully hops to 
a new minimum. The tapping strength A is then set to its 
initial value. The process is iterated until the rate of dis- 
covery of new global minima drops below some threshold 
value or a maximum number of iterations is reached. 

In the current implementation of the algorithm, the 
local minimization step is performed using the Fletcher- 
Reeves conjugate gradient algorithm*^. Analytic ex- 
pressions for the energy gradient and the Hessian ma- 
trix are coded in the program in order to reduce the 
number of evaluations of the objective function during 
the relaxation step to one single event. The Delaunay 
triangulation is calculated via the Dwyer's divide and 
conquer algorithm with alternate cuts^*, which runs in 
O(iVloglogiV) time, making the identification of the de- 
fects particularly fast. 

The main difference between TA and other hybrid algo- 
rithms (including Basin- Hopping) is that the escape move 
consists of adaptive displacements of defects only, rather 
than of the entire system. In the case of non-Euclidean 
crystals, where the conformation of the energy landscape 
is subtly related to the arrangement of topological de- 
fects, this mechanism is believed to explore the PES more 
accurately. In systems as Lennard- Jones clusters or spin- 
glasses, the PES is characterized by an exponential num- 
ber of local minima separated by energy barriers. For 
this reason the majority of the algorithms are specifi- 
cally designed to allow the system to overcome a barrier 
by providing it a significant amount of energy. If the 
energy landscape, however, is characterized by the pres- 
ence of multiple narrow funnels, as believed in this case, 
the previous methods become ineffective. A funnel rep- 
resents the basin of attraction of a given local minimum. 
If the global minimum is also located at the bottom of a 
funnel, an algorithm that is attempting to locate it via a 
sequence of local minimization steps has a chance to find 
it exclusively by starting from a configuration already at 
the muzzle of the funnel. Such possibility, however, is 
ruled out if all the particles are displaced simultaneously 
during the escape move and the system is abruptly moved 
to a completely different place in the energy landscape. 
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On the other hand, by adaptively tapping the defects it 
is possible to achieve a much finer inspection of the PES 



and possibly locate the funnel associated with the global 
minimum. A copy of our code is available by request. 
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